don <- read.table(file = "data/africa_pays_2018/data/africa_pays_2018.csv", # nom du fichier et chemin d'accès
sep = ";", # séparateur (ici, des points-virgule)
dec=",", # Type de décimale
header = TRUE, # ligne d'en-tête avec le nom des variables
encoding="UTF-8") # encodage adapté au françaisSTAT3 : REGRESSION
Comment modéliser une relation linéaire Y = a.X+b ?
Claude Grasland & Jean-Paul Nguesso
2025-06-01
1. PREPARATION DES DONNEES
1.1. Chargement du tableau principal
On charge le fichier des pays en 2018
1.2 Choix des deux variables à analyser
En dehors de iso3 et nom, on ne garde que deux variables que l’on renomme X et Y avec colnames() et que l’on convertit en type numérique général. Il suffira par la suite de modifier le choix des variables X et Y pour faire d’autres analyses.
CODE NOM X Y
1 AGO Angola 65.9 14.3
2 BDI Burundi 13.2 2.7
3 BEN Bénin 47.6 20.0
4 BFA Burkina Faso 29.7 16.0
5 BWA Botswana 69.8 47.0
6 CAF Rep. Centrafricaine 41.6 4.31.3 On est malin …
Mais comme on ne sait plus ce que sont X et Y, on le précise avec des chaînes de caractères qu’on pourra utiliser dans les graphiques. Et on peut préparer une version multilangue …
2. NORMALITE, VISUALISATION, CORRELATION
2.1 Vérification de la normalité de X et Y
La régression linéaire met en relation deux variables quantitatives X et Y dont on suppose que la distribution est normale (gaussienne) , c’est-à-dire unimodale et symérique.
Z<-rnorm(100000)
hist(Z, nclass=50, main="Loi normale (moyenne = 0 et écart-type = 1) ",border ="gray80",col="lightyellow",ylab=NULL,xlab=NULL,probability = TRUE,
cex.axis=0.9,ylim = c(0,0.55), cex.main=0.8,cex.lab=0.8,xlim=c(-3,3))
lines(density(Z,bw=0.2),col="red",lwd=1)
abline(v=c(0),col="blue",lty=1,lwd=1)
abline(v=c(-1,1),col="blue",lty=2,lwd=1)
abline(v=c(-2,2),col="blue",lty=1,lwd=1)
lines(c(-1,1),c(0.42,0.42))
text(0,0.45,"67%",cex=0.7)
lines(c(-2,2),c(0.5,0.5))
text(0,0.52,"95%",cex=0.7)2.1 Vérification de la normalité de X et Y
On peut tester la normalité des disributions par inspection visuelle à l’aide de hist()
On voit que X est assez symétrique et unimodale. En revanche Y est plutôt dissymétrique à gauche.
2.1 Vérification de la normalité de X et Y
Une analyse rapide avec boxplot() permet de repérer s’il y a des valeurs exceptionnelles.
2.1 Vérification de la normalité de X et Y
Une analyse rapide avec boxplot() permet de repérer s’il y a des valeurs exceptionnelles.
Pour X comme Y, on ne voit pas apparaître de valeurs exceptionnelles c’est-à-dire de points situés au delà de l’intervalle défini par la boîte à moustache.
2.1 Vérification de la normalité de X et Y
Les fonctions qqnorm() et qqline() sont plus précises …
X suit bien la droite et est proche d’une distribution gaussienne
2.1 Vérification de la normalité de X et Y
Les fonctions qqnorm() et qqline() sont plus précises …
Y ne suit pas bien la droite donc elle n’est pas tout à fait gaussienne.
2.1 Vérification de la normalité de X et Y
Mais la solution la plus précise est le test de Shapiro qui pose l’hypothèse H0 : la distribution est normale.
Shapiro-Wilk normality test
data: sel$X
W = 0.97546, p-value = 0.3931
Shapiro-Wilk normality test
data: sel$Y
W = 0.89916, p-value = 0.0005157Si le test est supérieur à 0.05, on peut considérer que la distribution est gaussienne ce qui est le cas de X. Dans le cas contraire (p <0.05) la distribution n’est pas gaussienne ce qui est le cas de Y.
Mais en pratique on accepte en général de faire des analyses de régression avec des distributions non-gaussienne, surtout s’il n’y a pas de valeurs exceptionnelles.
2.2 Visualisation de la forme de la relation
On peut faire un simple plot(X,Y). Mais on peut aussi créer pour cela une fonction personalisée adapté à ses préférences
monplot <- function (varX , varY, varN )
{
Ymin <- min(varY)
Ymax <- max(varY)
Ysup <- Ymax+(Ymax-Ymin)*0.1
plot(varX,varY,
main = titre, # titre
cex.main = 1, # police du titre
cex = 0.6, # taille des symboles
pch = 19, # cercles pleins
ylim = c(Ymin, Ysup),
col = "red") # couleur des symboles
text(varX,varY,varN,cex=0.5,pos=3) # nom des élément
abline(v=mean(varX),lty=2,lwd=1,col="blue") # moyenne X
abline(h=mean(varY),lty=2,lwd=1,col="blue") # moyenne Y
grid()
}2.2 Visualisation de la forme de la relation
Je peux désormais utiliser ma fonction monplot() !
2.2 Visualisation de la forme de la relation
Je peux décider de ne pas afficher le label des points.
2.3 Analyse de la corrélation
Je commence par calculer le coefficient de corrélation linéaire (r) et le pouvoir explicatif de X par rapport à Y (r2)
2.3 Analyse de la corrélation
Puis, je teste la significativité de la corrélation linéaire …
2.3 Analyse de la corrélation
… et je la compare à celle du coefficient de corrélation de rang de Spearman
Spearman's rank correlation rho
data: sel$X and sel$Y
S = 10474, p-value = 0.0007479
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.4656206 Les deux coefficients sont assez proches (+0.53 pour Pearson et +0.47 pour Spearman) ce qui est rassurant. En général, si les deux coefficients diffèrent il y a un problème de valeur exceptionnelle ou de relation non linéaire.
2.3 Analyse de la corrélation
On peut conclure des analyses précédentes que :
- il existe une relation significative (p-value < 0.01)
- cette relation est positive (r = +0.53 )
- cette relation a un pouvoir explicatif assez faible (r2 = 28%)
- la relation est monotone et globalement linéaire car le coefficient de Spearman est à peu près égale au coefficient de Pearson. Il n’y a donc pas besoin de transformer les variables ou retirer des valeurs exceptionnelles.
3. LA REGRESSION LINEAIRE
3.1 Hypothèses statistiques
Conditions a priori
- X et Y sont deux variables normales (gaussienne)
- il existe une corrélation significative entre X et Y (p< 0.05)
- X explique une part suffisamment forte de Y (r2 > 20% )
- Le nuage de point affiche une forme linéaire
- les points sont répartis de façon régulière le long du nuage de points
- Il n’y a pas de valeurs exceptionnelles susceptibles de perturber le calcul.
3.1 Hypothèses statistiques
Méthode des moindres carrés ordinaire (MCO)
La droite \(y_i = a.x_i + b + \epsilon_i\) qui minimise la somme des carrés des écarts entre les valeurs observées \(y_i\) et les valeurs estimées \(\hat{y_i}\) a pour équation :
\(COV(X,Y) = \sum_{i=1}^k \sum_{j=1}^k (x_{i}-\bar{x})^2.(y_{i}-\bar{y})^2\)
\(a = COV(X,Y) / (\sigma_X)^2\)
\(b = \bar{y} - a.\bar{x}\)
3.1 Hypothèses statistiques
Analyse de la variance
La somme des carré des écarts totale (\(SCE_{tot}\)) correspond à la variance de la variable à expliquer : \(SCE_{tot} = \sum_{i=1}^k (y_{i}-\bar{y})^2\)
La somme des carré des écarts résiduels (\(SCE_{err}\)) correspond à la somme des carrés des différences entre valeurs observées et estimées \(SCE_{error} = \sum_{i=1}^k (y_{i}-\hat{y})^2\)
Le pouvoir explicatif d’un modèle de régression correspond à la part de la variance de Y expliquée par X.
\(Var. expliquée = (SCE_{tot}-SCE_{res}) / SCE_{tot} = r(X,Y)^{2}\)
3.1 Hypothèses statistiques
Vérifications a posteriori
Un modèle de régression n’est valide que si les résidus de ce modèle \(\epsilon_i = (y_i - \hat{y}_i)\) remplissent les conditions suivantes :
- Normalité de la distribution des résidus
- Absence d’autocorrélation des résidus
- Homogénéité de la variance des résidus
- Absence de valeur à fort effet de levier
Si ces quatre conditions ne sont pas remplies, les estimations de Y en fonction de X seront entâchées d’erreur et leur intervalle de confiance ne sera pas valable.
3.2 La fonction lm()
La fonction lm() ou lm est l’abbréviation de linear model permet d’effectuer la plupart des modèles de régression linéaire basés sur la méthode des moindres carrés ordinaire. Sa syntaxe est a priori très simple et renvoie les coefficients b et a du modèle de régression.
3.2 La fonction lm()
Mais en réalité lm() crée une liste de résultats que l’on a intérêt à stocker pour en examiner les composantes une à une.
List of 12
$ coefficients : Named num [1:2] 1.122 0.534
..- attr(*, "names")= chr [1:2] "(Intercept)" "sel$X"
$ residuals : Named num [1:49] -22.042 -5.477 -6.562 -0.995 8.574 ...
..- attr(*, "names")= chr [1:49] "1" "2" "3" "4" ...
$ effects : Named num [1:49] -177.69 69.57 -3.8 2.06 10.97 ...
..- attr(*, "names")= chr [1:49] "(Intercept)" "sel$X" "" "" ...
$ rank : int 2
$ fitted.values: Named num [1:49] 36.34 8.18 26.56 17 38.43 ...
..- attr(*, "names")= chr [1:49] "1" "2" "3" "4" ...
$ assign : int [1:2] 0 1
$ qr :List of 5
..$ qr : num [1:49, 1:2] -7 0.143 0.143 0.143 0.143 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:49] "1" "2" "3" "4" ...
.. .. ..$ : chr [1:2] "(Intercept)" "sel$X"
.. ..- attr(*, "assign")= int [1:2] 0 1
..$ qraux: num [1:2] 1.14 1.27
..$ pivot: int [1:2] 1 2
..$ tol : num 1e-07
..$ rank : int 2
..- attr(*, "class")= chr "qr"
$ df.residual : int 47
$ xlevels : Named list()
$ call : language lm(formula = sel$Y ~ sel$X)
$ terms :Classes 'terms', 'formula' language sel$Y ~ sel$X
.. ..- attr(*, "variables")= language list(sel$Y, sel$X)
.. ..- attr(*, "factors")= int [1:2, 1] 0 1
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:2] "sel$Y" "sel$X"
.. .. .. ..$ : chr "sel$X"
.. ..- attr(*, "term.labels")= chr "sel$X"
.. ..- attr(*, "order")= int 1
.. ..- attr(*, "intercept")= int 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
.. ..- attr(*, "predvars")= language list(sel$Y, sel$X)
.. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
.. .. ..- attr(*, "names")= chr [1:2] "sel$Y" "sel$X"
$ model :'data.frame': 49 obs. of 2 variables:
..$ sel$Y: num [1:49] 14.3 2.7 20 16 47 4.3 46.8 23.2 8.6 8.7 ...
..$ sel$X: num [1:49] 65.9 13.2 47.6 29.7 69.8 41.6 51 56.7 44.8 67.2 ...
..- attr(*, "terms")=Classes 'terms', 'formula' language sel$Y ~ sel$X
.. .. ..- attr(*, "variables")= language list(sel$Y, sel$X)
.. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
.. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. ..$ : chr [1:2] "sel$Y" "sel$X"
.. .. .. .. ..$ : chr "sel$X"
.. .. ..- attr(*, "term.labels")= chr "sel$X"
.. .. ..- attr(*, "order")= int 1
.. .. ..- attr(*, "intercept")= int 1
.. .. ..- attr(*, "response")= int 1
.. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
.. .. ..- attr(*, "predvars")= language list(sel$Y, sel$X)
.. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
.. .. .. ..- attr(*, "names")= chr [1:2] "sel$Y" "sel$X"
- attr(*, "class")= chr "lm"3.3 Analyse du modèle
Un résumé des résultats principaux est fourni avec summary() appliqué à l’objet créé par lm().
Call:
lm(formula = sel$Y ~ sel$X)
Residuals:
Min 1Q Median 3Q Max
-28.337 -11.417 -2.522 12.700 33.105
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.1223 6.1037 0.184 0.855
sel$X 0.5344 0.1244 4.295 8.68e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 16.2 on 47 degrees of freedom
Multiple R-squared: 0.2819, Adjusted R-squared: 0.2666
F-statistic: 18.45 on 1 and 47 DF, p-value: 8.68e-05On obtient ainsi :
- l’équation de la droite Y = a.X+b
- la significativité et l’intervalle de confiance de a et b
- le pouvoir explicatif du modèle \(r(X,Y)^2\)
3.4 Tracé de la droite de régression
On peut tracer la droite de régression avec abline()
Cette droite a pour équation Y = 0.5X + 2.8 ce qui veut dire que chaque fois que l’urbanisation augemente de 1 pt l’usage d’internet augmente de 0.5 pts. Un pays qui aurait une urbanisation de 50% devrait donc avoir un taux d’usage d’internet égal à 0.5 x 50 + 2.8 = 27.8%
3.5 Analyse des résidus
On peut extraire de l’objet créé par lm() les valeurs estimées de Y et les résidus c’est-à-dire les erreurs d’estimation.
sel$Y_estim<-monmodel$fitted.values # Valeurs estimées
sel$Y_resid<-monmodel$residuals # Valeurs résiduelles
sel <- sel[order(sel$Y_resid),] # Tri du tableau en fonction des résidus
sel CODE NOM X Y Y_estim Y_resid
10 COG Congo 67.2 8.7 37.036626 -28.3366260
40 SOM Somalie 45.3 2.0 25.332411 -23.3324109
24 LBY Libye 80.3 21.8 44.037777 -22.2377775
1 AGO Angola 65.9 14.3 36.341855 -22.0418552
14 ERI Erythrée 40.4 1.3 22.713660 -21.4136596
23 LBR Libéria 51.4 8.0 28.592489 -20.5924891
20 GNB Guinée-Bissau 43.6 3.9 24.423865 -20.5238646
6 CAF Rep. Centrafricaine 41.6 4.3 23.354986 -19.0549865
9 COD Congo, Rép. dém. du 44.8 8.6 25.065191 -16.4651914
39 SLE Sierra Leone 42.3 9.0 23.729094 -14.7290938
19 GMB Gambie 61.6 19.8 34.043767 -14.2437674
21 GNQ Guinée équatoriale 72.4 26.2 39.815709 -13.6157090
27 MDG Madagascar 37.6 9.8 21.217230 -11.4172303
43 TGO Togo 42.0 12.4 23.568762 -11.1687621
28 MLI Mali 42.8 13.0 23.996313 -10.9963133
29 MOZ Mozambique 36.3 10.0 20.522460 -10.5224596
48 ZMB Zambie 43.8 14.3 24.530752 -10.2307524
30 MRT Mauritanie 54.1 20.8 30.035475 -9.2354745
8 CMR Cameroun 56.7 23.2 31.425016 -8.2250160
42 TCD Tchad 23.2 6.5 13.521308 -7.0213081
3 BEN Bénin 47.6 20.0 26.561621 -6.5616207
2 BDI Burundi 13.2 2.7 8.176918 -5.4769176
33 NER Niger 16.5 5.3 9.940566 -4.6405665
37 SSD Sud Soudan 19.8 8.0 11.704215 -3.7042153
18 GIN Guinée 36.3 18.0 20.522460 -2.5224596
4 BFA Burkina Faso 29.7 16.0 16.995162 -0.9951619
22 KEN Kenya 27.3 17.8 15.712508 2.0874918
31 MWI Malawi 17.1 13.8 10.261230 3.5387701
45 TZA Tanzanie 34.2 25.0 19.400138 5.5998624
15 ETH Ethiopie 21.0 18.6 12.345542 6.2544578
17 GHA Ghana 56.4 39.0 31.264684 7.7353157
5 BWA Botswana 69.8 47.0 38.426168 8.5738325
49 ZWE Zimbabwe 32.2 27.1 18.331259 8.7687405
46 UGA Ouganda 24.1 23.7 14.002303 9.6976968
36 SDN Soudan 34.8 30.9 19.720801 11.1791990
35 RWA Rwanda 17.3 21.8 10.368118 11.4318823
25 LSO Lesotho 28.4 29.0 16.300391 12.6996089
11 DJI Djibouti 77.9 55.7 42.755124 12.9448762
16 GAB Gabon 89.6 62.0 49.008061 12.9919394
34 NGA Nigéria 50.8 42.0 28.271826 13.7281743
7 CIV Côte d'Ivoire 51.0 46.8 28.378713 18.4212865
47 ZAF Afrique du Sud 66.7 56.2 36.769406 19.4305935
38 SEN Sénégal 47.5 46.0 26.508177 19.4918232
12 DZA Algérie 72.9 59.6 40.082929 19.5170714
32 NAM Namibie 50.5 51.0 28.111494 22.8885060
13 EGY Egypte 42.7 46.9 23.942869 22.9571306
44 TUN Tunisie 69.1 64.2 38.052060 26.1479398
26 MAR Maroc 62.8 64.8 34.685094 30.1149058
41 SWZ Swaziland 23.9 47.0 13.895415 33.10458463.5 Analyse des résidus
Les résidus négatifs correspondent aux pays qui ont moins accès à internet que ce que leur urbanisation laisserait prévoir :
CODE NOM X Y Y_estim Y_resid
10 COG Congo 67.2 8.7 37.03663 -28.33663
40 SOM Somalie 45.3 2.0 25.33241 -23.33241
24 LBY Libye 80.3 21.8 44.03778 -22.23778
1 AGO Angola 65.9 14.3 36.34186 -22.04186
14 ERI Erythrée 40.4 1.3 22.71366 -21.41366- Ainsi le Congo n’a que 8.7% d’accès à internet (Y) alors que son taux d’urbanisation est de 67.2% (X) ce qui laisserait prévoir un taux d’accès à internet de 37.0% (Y_estim). Il a donc un résidus de -28 pts (Y_resid) qui signifie que son taux d’accès à internet est beaucoup plus faible que ce que laisserait prévoir son urbanisation.
3.5 Analyse des résidus
Les résidus positifs correspondent aux pays qui ont plus accès à internet que ce que leur urbanisation laisserait prévoir :
CODE NOM X Y Y_estim Y_resid
32 NAM Namibie 50.5 51.0 28.11149 22.88851
13 EGY Egypte 42.7 46.9 23.94287 22.95713
44 TUN Tunisie 69.1 64.2 38.05206 26.14794
26 MAR Maroc 62.8 64.8 34.68509 30.11491
41 SWZ Swaziland 23.9 47.0 13.89542 33.10458-Ainsi le Swaziland a 47 % d’accès à internet (Y) alors que son taux d’urbanisation est de 23.9% (X) ce qui laisserait prévoir un taux d’accès à internet de 13.9% (Y_estim). Il a donc un résidus de +33 pts (Y_resid) qui signifie que son taux d’accès à internet est beaucoup plus fort que ce que laisserait prévoir son urbanisation.
4. COMPLEMENTS STATISTIQUES (facultatifs)
4.1 Hypothèses statistiques
Vérifications a posteriori
Un modèle de régression n’est valide que si les résidus de ce modèle \(\epsilon_i = (y_i - \hat{y}_i)\) remplissent les conditions suivantes :
- Normalité de la distribution des résidus
- Absence d’autocorrélation des résidus
- Homogénéité de la variance des résidus
- Absence de valeur à fort effet de levier
Si ces quatre conditions ne sont pas remplies, les estimations de Y en fonction de X seront entâchées d’erreur et leur intervalle de confiance ne sera pas valable.
On charge le package car (companion to applied regession) pour réaliser certains des diagnostics
4.1 Diagnostic 1 : Indépendance des résidus ?
L’objectif est de savoir si les résidus se répartissent régulièrement de part et d’autre de la droite de régression tout au long de celle-ci. Si c’est bien le cas le graphique residuals Vs Fitted généré par plot(monmodel,1) devrait donner une droite horizontale :
4.1 Diagnostic 1 : Indépendance des résidus ?
4.1 Diagnostic 1 : Indépendance des résidus ?
On peut tester statistiquement l’indépendance des résidus à l’aide du test de Durbin-Watson qui mesure si deux valeurs successives ont des résidus proches. L’indépendance des résidus est rejetée si p-value < 0.05
lag Autocorrelation D-W Statistic p-value
1 0.00312332 1.94811 0.856
Alternative hypothesis: rho != 0Ici on trouve p-value > 0.05 donc les résidus sont indépendants.
Diagnostic 2 : Normalité des résidus ?
L’objectif est de savoir si les résidus ont une distribution normale Si c’est bien le cas le graphique généré par plot(monmodel,2) devrait donner une droite oblique :
Diagnostic 2 : Normalité des résidus ?
Diagnostic 2 : Normalité des résidus ?
On peut tester statistiquement la normalité des résidus à l’aide du test de Shapiro-Wilk. Les résidus sont normaux si p-value > 0.05
Shapiro-Wilk normality test
data: monmodel$residuals
W = 0.96515, p-value = 0.1543Ici on trouve une p-value supérieure à 0.05 donc la distribution des résidus est approximativement gaussienne.
Diagnostic 3 : Homogénéité des résidus ?
L’objectif est de savoir si la variance des résidus est constante, c’est-à-dire si il s’écarte environ de la même distance tout au long de la droite . Si c’est bien le cas le graphique généré par plot(monmodel,3) devrait donner une droite horizontale
Diagnostic 3 : Homogénéité des résidus ?
Diagnostic 3 : Homogénéité des résidus ?
On peut tester statistiquement l’homogénéité des résidus à l’aide du test de Breush-Pagan. L’hypothèse d’homogénéité est rejetée si la p-value est inférieure à 0.05.
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 3.029338, Df = 1, p = 0.081771Ici, la p-value est supérieure à 0.05 donc les résidus sont globalement homogènes.
Diagnostic 4 : Absence de valeur exceptionnelles ?
L’objectif est de savoir s’il existe des valeurs qui exercent une influence exceptionnelle sur les résultats de la régression. On peut reprérer ces valeurs de plusieurs manières, notamment à l’aide de la distance de Cook générée par plot(monmodel,4).
Diagnostic 4 : Absence de valeur exceptionnelles ?
On repère les pays les plus susceptibles de perturber le modèle (Sierra Léone, Sud Soudan et Côte d’Ivoire)
Diagnostic 4 : Absence de valeur exceptionnelles ?
Le test statistique de Bonferroni permet de déterminer s’il existe des valeurs exceptionnelles avec une p-value < 0.05.
No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:
rstudent unadjusted p-value Bonferroni p
CIV 2.176162 0.034715 NAIci, on voit que la Côte d’Ivoire est le seul pays ayant une influence sur le modèle. Mais la p-value est tout juste inférieure à 0.05 donc ce n’est pas très grave.